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We discuss three formally different formulas for normalization of quasinormal modes currently in 
use for modeling optical cavities and plasmonic resonators and show that they are complementary 
and provide the same result. Regardless of the formula used for normalization, one can use the 
norm to define an effective mode volume for use in Purcell factor calculations. 


I. INTRODUCTION 


Optical cavities and plasmonic resonators enable tight 
localization of electromagnetic Helds to length scales on 
the order of the wavelength or shorter. In addition, the 
resonant nature of these material systems suggest a dis¬ 
crete nature of the possible electromagnetic waves, sim¬ 
ilar to the bound states of the hydrogen atom. The 
analogy is somewhat deceptive because of the inherently 
leaky nature of electromagnetic resonators, which leads 
to a continuous dissipation of energy via radiation or 
absorption in the material. Nevertheless, distinct reso¬ 
nances show up as peaks in spectra from scattering calcu¬ 
lations or measurements, for example, in which the width 
of the peak is related to a finite lifetime due to the dissipa¬ 
tion of energy; the associated scattered fields often show 
signatures of an underlying modal structure. Throughout 
the literature, these resonances are frequently referred to 
as “modes”, but an explicit mathematical definition of 
the modes is rarely given, and this may cause both con¬ 
ceptual and practical difficulties when using the modes 
in calculations. 
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FIG. 1: (Color online) Left: Illustration of a gold nanorod 
dimer. Right: Absolute value of the fundamental QNM f c (r) 
in the plane z = 0. The QNM is scaled to unity at the position 
r c = (0, 0, 0) in the center between the two rods. 


In computational nanophotonics, the optical or plas¬ 
monic response of resonant systems are often modeled by 
use of time-domain calculations with Perfectly Matched 
Layers (PMLs) and a run-time Fourier transform to ex¬ 
tract the modal field distribution. This technique has 
been very successful for modeling high quality optical 
cavities for which the spectral resonances are well sep¬ 
arated. In particular, the modes so calculated have 
been used for estimates of the spontaneous emission rate 
enhancement via the Purcell formula 1 } from which the 
spontaneous emission rate is expected to scale with the 
ratio of the cavity quality factor Q to the effective mode 
volume V e g. The effective mode volume describes the 
degree of light localization in the cavity and has a clear 
definition for closed resonators in terms of the integrated 
energy density of light in the mode. The definition, how¬ 
ever, cannot be directly extended to dissipative systems 
for which the mode extends over all space, and this has 
conceptual consequences for application of the Purcell 
formula which manifestly pertains to dissipative systems 
with a finite Q value. Theoretically, the response of dissi¬ 
pative resonant structures can be understood in terms of 
so-called quasinormal modes 2 ® (QNMs), which appear 
as discrete solutions to the wave equation and satisfy a 
radiation condition compatible with the propagation of 
light away from the resonator. Figure [T] shows an exam¬ 
ple of a plasmonic resonator, in the form of a nanorod 
dimer, and the field profile of the fundamental QNM in 
the dimer. The use of PMLs is one possible realization 
of a radiation condition for resonators in homogeneous 
media, wherefore it was pointed out in Ref.[5]that the op¬ 
tical cavity modes commonly calculated by time-domain 
techniques are exactly the QNMs, and that the Purcell 
formula can be derived within a QNM framework such as 
developed in Refs.lSHIUl for example; this leads naturally 
to a definition of the effective mode volume which is com¬ 
patible with the leaky nature of the modes. The effective 
mode volume is intimately connected with the normal¬ 
ization of the QNMs, and two alternative formulations 
of the effective mode volume were derived in Refs. EH 
and [H respectively, based on different formulations of 
the norm. In this article we show how the three different 
normalization integrals are closely related and provide 
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the same result. In essence, we show that the differ¬ 
ences in the formulas can be understood as alternative 
regularization procedures for an inherently ill-behaved 
integral. Although we illustrate the formal connection 
between the three formulations, it is worth emphasizing 
that they were derived by use of very different approaches 
in Refs. is no and ESI The equivalence of the formulas, 
therefore, should be regarded as a profound strength of 
the entire modeling framework for localized electromag¬ 
netic resonators based on QNMs. 

The article is organized as follows: In section [Tl] we 
introduce the three formulations of the QNM norm and 
discuss how they relate to each other. Section [Til] pro¬ 
vides explicit example calculations to illustrate the prac¬ 
tical evaluation as well as the equivalence of the norms. 
Section m summarizes and concludes the article. 


II. NORMALIZATION INTEGRALS 

We consider electric field QNMs f M (r) of resonators 
embedded in a homogeneous background defined as the 
eigenmodes of Maxwell’s wave equation 

VxVxf^rJ-^uJf^^O, (1) 

in which = uj^/c is the ratio of the angular resonance 
frequency to the speed of light and e r (r, <D M ) is the rela¬ 
tive permittivity. Because of the leaky nature of the res¬ 
onators, the differential equation should be augmented 
by a radiation condition to only allow outgoing wave so¬ 
lutions at large distances. For the case of a homogeneous 
environment, we argue that this radiation condition is 
formally the Silver-Miiller radiation condition^ in the 
form 

r - - - 

- x V x f^(r) -a -m B ^f^( r ) as r — >■ oo, (2) 

where ne is the refractive index of the homogeneous back¬ 
ground material. The Silver-Miiller radiation condition 
is commonly imposed to correctly model scattered fields. 
Therefore, we argue that this is also the radiation con¬ 
dition for the QNMs, which can be considered solutions 
to the scattering problem with no source. Because it is 
defined only in the limit of infinite distances, Eq. <[2j) is 
rarely useful in numerical calculations of QNMs. Instead, 
one can use PMLs to mimic the radiation condition by 
removing reflections from the calculation domain bound¬ 
ary, or one can use calculation methods based on the 
Green tensor which manifestly respects the radiation con¬ 
dition. As a result of the radiation condition, the differ¬ 
ential equation problem becomes non-Hermitian, where¬ 
fore a number of well-known text book results for Hermi- 
tian eigenvalue problems do not apply. In particular, the 
discrete resonance frequencies - iy M are com¬ 

plex with a negative imaginary part, from which the Q 
value can be calculated as Q^ The complex 

resonance frequencies in combination with the radiation 
condition leads to a divergent behavior of the QNMs as 


a function of distance from the resonator and makes the 
normalization non-trivial. As discussed in the introduc¬ 
tion, the problem of normalization for QNMs has been 
addressed by at least three different methods leading to 
three different formulas that we discuss in detail below. 
Before moving on, however, it is worth elaborating a bit 
on the connection between the mathematical concept of 
the norm, and the more physical concept of mode volume. 


Regardless of the formula used for normalization, one 
can use the norm to define a generalized effective mode 
volume for the QNM of interest (denoted by n = c) a^ 


((fjfc)) 

e r (r c )f c 2 ( r c)’ 


( 3 ) 


where ((f c |f c )} denotes the QNM norm (in any formula¬ 
tion compatible with the divergent behavior of the QNMs 
at large distances) and f(?(r c ) = f c (r c ) • f c (r c ). The po¬ 
sition r c denotes a point of special interest to the given 
resonator. For optical cavities, this will often be the cen¬ 
ter of the cavity in which the field is largest, and the 
generalized effective mode volume then gives a measure 
for the maximum light-matter interaction that one can 
obtain in the cavity. For general resonators there may 
be no unique way of defining the center, or it may be 
that the mode vanishes at this point. For plasmonic res¬ 
onators, in particular, the field may be several orders of 
magnitude larger at the edges than within the material. 
Despite the large field values, however, the edges are of¬ 
ten uninteresting for light-matter interaction because of 
enormous non-radiative decay and quenching. In general, 
therefore, we specify r c explicitly for each resonator. The 
generalized effective mode volume is complex in general, 
and is related to the effective mode volume V e g of the 
Purcell formula as 


1 



( 4 ) 


This definition follows directly and quite naturally from 
a derivation of the original Purcell formula in the frame¬ 
work of QNMs, either using a Green tensor approacll^l- R ^ I 
or a mode expansion formulation 1 ^. It follows from these 
derivations, that the usefulness of the effective mode vol¬ 
ume for Purcell factor calculations is restricted to mate¬ 
rial systems or frequency ranges for which a single mode 
dominates the response. Moreover, the use of Eq. Q in 
Purcell factor calculations is valid in the limit w c /u; c ~ 1 
which may not hold for very low-Q optical cavities or 
plasmonic resonators. In cases where it does not hold, 
one can still recover the Purcell formula by redefining 
either Eq. ([3]) or Q, but in practice it is often easier 
and more convenient to calculate the spontaneous emis¬ 
sion enhancement directly from a QNM expansion of the 
Green tensor. For plasmonic structures, the usefulness 
of the original Purcell formula is further limited by non- 
radiativ e dec ay channels at positions very close to the 
material— .11 ^1. 
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A. Normalization by Lai et al. 


Lai et al. introduced a normalization for QNMs, which 
was applied to the study of spherically symmetric and 
non-dispersive material systems in Ref. Rewriting 
slightly, we express this norm as 

((f M |f/x))Lai = Jim { [ e r (r)f M (r) • f M (r)dR 

^->■00 l Jy 

+ [ *U r ) ' ^0)^}, ( 5 ) 

JdV J 

where V is a spherical volume with boundary dV and dR 
and dA are differential volume and area elements, respec¬ 
tively. The normalization in Ref. [S] appears to date back 
to early work by Zel’dovich on the theory of unstable 
states 1 A It was later adopted, in the form of Eq. ©I, to 
introduce a generalized effective mode volume for leaky 
optical cavitieiP. The formula can be extended to disper¬ 
sive resonators in homogeneous and non-dispersive sur¬ 
roundings by the substitution 7 e r (r) -A <r(r,u; M ), where 

cr( r , w) =* (w 2 e r (r, w)). (6) 


The norm in Eq. © suffers from a pathological defi¬ 
ciency in the case of degenerate modes of high-symmetry 
material systems. For spherical resonators, for example, 
one can choose the azimuthal dependence of the QNMs 
to be of the form exp{±irw^}, for integer values of m. 
This choice will lead to a vanishing norm. As discussed 
in Ref. (6j it may therefore be more appropriate to for¬ 
mulate the norm in terms of the adjoint fields which, 
for spherical systems, is obtained by complex conjuga¬ 
tion of the angular dependence onljP It is clear, how¬ 
ever, that any perturbative breaking of the symmetry will 
lift the degeneracy and lead to well-defined QNMs for 
which the norm applies immediately. Therefore, in such 
cases where the norm vanishes, one can always choose 
linear combinations of the degenerate QNMs for which 
the norm applies as written in Eq. ©■ In the case of 
the spherical resonator, for example, one can choose the 
azimuthal dependencies of the QNMs to be of the form 
exp{im<p} ± exp{— imip}; in this case Eq. ([5]) applies di¬ 
rectly. 


For any resonator, we can evaluate the norm in Eq. © 
by first integrating over a spherical domain Vg of ra¬ 
dius Rq completely enclosing the resonator; this integral 
is well behaved and can be easily evaluated. Next, for 
the region r > Rg, we can write the QNMs in terms of 
spherical Hankel functions hiin^k^r) and vector spheri¬ 
cal harmonicf®. Therefore, the norm may be rewritten 
in the form 


««>Lai = / e r (r)f M (r).f M (r)dR + ^if/JL (7) 

J Vn i 

u Im 


in which I[ = lim r^ > . 00 {I[(R)} with 

,-R 


I[(R) = I n‘^hf(nBk f _ l r)r 2 dr + i^-hf(nBk f j,R)R 2 , 

2 


' Ro 


( 8 ) 


and l(* n includes the angular integration of the vector 
spherical harmonics, which is independent of R (and fi¬ 
nite). Here, and in the following, we use linear combi¬ 
nations of vector spherical harmonics of the form P^ n = 
Pirn ± P* m , where P i m is any vector spherical harmonic 
of order ( l,m ) in the formulation of Ref. ITS} The ar¬ 
gument, originally by Lai et alW. is now, that because 
of the outgoing wave nature of the Hankel functions at 
large distances, the increase in the volume integral is ex¬ 
actly balanced by the additional surface integral. For 
large arguments, the spherical Hankel functions tend to 
the limiting form 


h t (z) -Ae- i7r(z+1)/2 -e iz , 


(9) 


so that inserting in Eq. 
spect to R we find 

1 


and differentiating with re- 


d R I[{R) = (-1)('+!) + — 

fc 2 1 2nsfcju v 2 J 

= 0 , ( 10 ) 


suggesting that one can assign a well-defined value to the 
norm as R —> oo. 


In practice, direct application of Eq. ([5]) leads to an 
integral that seems to quickly converge towards a finite 
value, but in fact oscillates about this value with an 
amplitude that eventually starts to grow (exponentially) 
with the distance, albeit slowly compared to the length 
scales in typical calculations. This was noted in Ref. 3 
where the oscillations were observed only for the cavity 
with the lowest quality factor (Q « 16). The source of 
the oscillating integral can be traced back to the complex 
resonance frequency, which means that the wavenumber 
kfj, is complex as well. Evaluation of Eq. (|8| then leads 
to residual radius dependent terms of the form 

™ = ( 11 ) 


where P n (R) denotes general polynomials in R of order 
n > 1; the surface term in Eq. © cancels the first or¬ 
der term. Thus, while Eqs. © and (10) appear to be 
formally correct also for complex arguments, the limit 
R —> oo in practice leads to a position dependent phase 
difference between the Hankel function and its limiting 
form, which makes the limit non-trivial to perform along 
the real axis. In principle, however, it is possible to reg¬ 
ularize the integral using coordinate transforms; a tech¬ 
nique that has also been also for normalization of leaky 
modes in waveguides 19 . Formally, we can rewrite the in- 
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tegral as 


as 


U rR 

nfthf(nBkur)r 2 dr + \ Tl ^-h 2 {nBknR)R 2 

Ro 2 

- [ R d r fS i (r)dr+ /* dJl^r) dr} 

J-Ro J Ro J 

= Fi(R 0 ) + f^Ro) + lim f d r flZ(r)dr, (12) 

JRn 


where 

-F}( r ) =-2^ (hU n Bk fJ .r) - hi- 1 (n B k ll r)hi +1 (nBk IJ ,r)) 

(13) 

is minus the antiderivative of /z(r) = {n-e,k^r)r 2 . In 

the last integral, we can now change the integration con¬ 
tour by writing r = i ? 0 + to find If = Fi(R 0 ). For the 
case of a spherical resonator, we provide an example of 
this approach in section III A Obviously, this regulariza¬ 
tion approach is rather cumbersome unless the expansion 
in spherical wave functions is explicitly known; also, it is 
rarely necessary in practice. Although the integral oscil¬ 
lates as a function of calculation domain size, the asso¬ 
ciated uncertainty in the calculated norm is much lower 
than the general numerical uncertainty for most QNM 
calculations. We elaborate on this issue in sections IlII Bl 
and lHICl 


B. Normalization by Sauvan et al. 

Sauvan et a/P-^ used a different approach based on the 
Lorentz reciprocity theorem to introduce a normalization 
for QNMs, in general dispersive and possibly magnetic 
materials, that can be written as 


«W»Sa 


1 


1 /i| J -/i//bauvan o / M 
1 JV 


f M (r) 

- g M (r) • At(r,w)g M (r)dV, (14) 

where g M (r) is the magnetic field QNM satisfying V x 
f)i(r) = iw^oMr(i\w M )g M (r), 77 (r,w) = 3 w (we r (w,r)) and 
k(t ,w) = d u {uiyLQ^ T {uj 1 r))/eo; eo and /z 0 denote the per¬ 
mittivity and permeability of free space, respectively, and 
/x r (r, ui) is the (possibly dispersive) relative permeability. 
This normalization seems to originate from the theory of 
leaky modes in waveguide^®. In the original formulation 
and implementation in Ref. 111) there is no explicit re¬ 
quirement of the volume to tend to infinity. Instead, the 
integral is performed over both the usual calculation do¬ 
main and the surrounding PMLs in which the real space 
coordinates are rotated into the complex plane, causing 
the propagating waves to decrease exponentially. 

For isotropic and non-magnetic materials, we may 
write Eq. (14) in terms of the electric field QNMs only 


«^|f M ))sauvan = 77 J f?(r, w^f^r) • f M (r) 

+ i(Vxf M (r)).(Vxf,(r))dV. (15) 

V 

For such materials, the correspondence with the normal¬ 
ization of Lai et al. was shown in Ref. m and is repeated 
here for completeness. First, we use the vector general¬ 
ization of Green’s identity of the first kincPS, 

/ (V x P) • (V x Q) - P • V x V x Q dV 


/ n • (P x V x Q) d.A, (16) 

JdV 


to rewrite Eq. (15) as 


((f/.|f M ))sauvan J v( r,^)f^(r)'f„(r) 


+ J^(r)-Vx Vxf M (r)dF 
^[1 

+ ^/ w n -( f " <r,xVxf - (r) ) dA 

(17) 

Then, using the wave equation and Eq. ([ 2 ]) in the limit 
V —> 00, we recover Eq. ([ 5 ]) with the substitution e r (r) —► 
cr(r, tip). Importantly, the integrand in Eq. ( |14| is invari¬ 
ant under the coordinate transformations of the PMLsP. 
Therefore, the integral must have a well-defined value 
also under the trivial transformation where n o c oordi- 
nate rotation is performed^. In view of Eq. ( fl7| >, this 


argument is a complement to Eq. (101 in suggesting that 
the norm has a well-defined value as R —» 00. 


C. Normalization by Muljarov et al. 


Muljarov et a/P^ 13 derived a normalization which is 
similar in structure to Eq. ([ 5 ]), but with a different surface 
term, 

((f/df/^Muijarov = / cr(r, w M )f M (r) • f)dr)dr 
Jv 

+ i/ W r ) ' 9 s (r9 r f M (r)) 

JdV 

— r(d r f M (r)) ■ (d s f M (r))dA, (18) 

where er(r,u;) is dehned in Eq. © , and d s denotes the 
derivative in the direction normal to the surface of the 
volume V, which needs not be spherical and needs not 
extend to infinity. Comparing to the two other formulas, 
Eq. (181 makes explicit reference to a spherical coordinate 


system, the center of which one is left free to choose. 

For the case of a spherical domain of radius Rq, we 
can show the correspondence with the normalization of 
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Lai et al. by first evaluating the integrals I[ in Eq. (12) 
using the recurrence relations 


hi-i(kr) = l -^^hi(kr) + \d r hi(kr) (19) 

hr k 

hi+i(kr) = -khi(kr) - ^ d r hi(kr ), (20) 

and the defining equation for the spherical hankel func¬ 
tions, r 2 d rr hi(r) + 2 rd r hi(r) + (r 2 — l(l +1 ))hi(r) = 0, to 
find 

2 

Fi(r) = (^hi{nBk^,r)d r (rd r hi {n^k^r)) 

- r(9 r /i z (n B fc At r )) 2 )- ( 21 ) 


Inserting in Eq. (|7|) and rearranging the terms by mul¬ 
tiplying onto the spherical harmonics, we can write the 
angular integration as a surface integral in terms of the 
QNMs in the exact form of Eq. (18). 


III. EXAMPLE CALCULATIONS 


In this section we present a number of example calcu¬ 
lations to illustrate that the three formulations give the 
same value for the norm in practical calculations. First, 
we consider the case of a homogeneous sphere, for which 
the functional form of the QNMs are known analytically 
and we can carry out the regularization of Eq. © explic¬ 
itly. This example was also considered by Muljarov et al. 
in Ref. Qj], Next, we revisit the photonic crystal cavity 
of Ref. E| and the plasmonic nanorod dimer of Ref. 1221 to 
elaborate on the practical evaluation of Eq. ([5]) . 


A. Homogeneous sphere 


As a first set of example applications of the three for¬ 
mulations, we consider QNMs of homogeneous spheres 
with radius a in air. Following Muljarov et al., we con¬ 
sider electric field QNMs of the transverse kind which 
may be written, for positions outside the sphere, in terms 
of the vector spherical harmonics 3>z m (0, ip) as 

Mr) = MV-)-^(d± (22) 


with m > 0, where ^i m {0,p) = r x VYi m (0,ip), and 
Yim(0,<p) denote the (scalar) spherical harmonics. For 
simplicity, we limit the discussion to the case of purely 
real angular dependence (corresponding to the “+” in 
Eq. (22)), and 1 = 1 for which the spherical Hankel func¬ 
tion may be written explicitly as 


hi{z) 


-(* + i) 


(23) 


Inserting in Eq. ([5]) it follows that 

«UU> M - ns--£)«■*• 

+ lim {flZm , (24) 


where 


= r [ u(r)Wr) • fim(r)dV (25) 
Jo Jn 

denotes the integral over the finite sphere, and the resid¬ 
ual R dependent factor is 




i c 2ik u R 

j&R 2 

h 1 


(26) 


By rewriting the norm as in Eq. (12) we can regularize 
the integral and conclude that the norm is 


«fi m |fim»Lai = /[<“ + 2 (4- - J- | e 2i ^“. (27) 

\2fc3 k%a) 

In a similar way, we can evaluate Eq. © to find that 


((flm|flm))s. 


_ T r < a O 
auvan — ± l m ' ^ 


~2i k„. 


2kl kf l a / 


+ Jim {/w M), 

it—>-00 


where 


/Sauvan(^) = 


2i 


J2\k u R 


kfR 3 kf,R 2 

\ r 1 ' / 


(28) 


(29) 


Clearly, Eq. (28) can be regularized in a way completely 
analogous to Eq. (241, which is nothing but a different 
way of introducing the coordinate transform that was 
performed by use of PMLs in Ref. [Tl] With this ap¬ 
proach, one finds immediately that the two norms are 
identical. 

Finally, using Eq. (181 one can directly verify that 


((flm|flm))Muljarov — 4V + 2 


2fc3 k*a y 


^2i k u a 


(30) 


so that it equals the result of the other two formulations 
of the norm — notably, without an explicit need for reg¬ 
ularization. 


B. Photonic crystal cavity 

As a second example, we consider the simple two di¬ 
mensional photonic crystal cavity of Ref. [5> which is 
formed by six cylinders of relative permittivity e rot is = 
11.4 in air. The cylinders are arranged in a hexagon 
with side length a and have radii r = 0.15 a. We con¬ 
sider the out-of-plane polarization, for which the cav- 
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ity supports a QNM f c (r) with complex resonance fre¬ 
quency uj c cl/2ttc = 0.425862 — 0.013539i, corresponding 
to Q ~ 16. Figure [5] shows the real part of f c (r) which is 
scaled to unity in the cavity center r c = (0, 0). 



Position x/a 


FIG. 2: (Color online) Real part of the fundamental QNM 
f c (r) in a small cavity made from high-index rods arranged 
in a hexagon with side length a, as indicated by black circles. 
The QNM is scaled to unity in the cavity center r c = (0,0). 

For these QNM calculations, we used a Fredholm type 
integral equation in which the QNMs appear as the self- 
consistent solutions to a scattering problem with no in¬ 
cident field 5 . The equation was solved by an iterative 
procedure in which, for each guess of a resonance fre¬ 
quency Wg U e SS , we set up the integral equation and solved 
it to find the eigenvalue closest to Wg U e SS ; the iterations 
then continued until the difference was smaller than some 
prescribed tolerance. Specifically, for each w guess , we 
discretized the integral eigenvalue equation by expand¬ 
ing the background Green tensor and the unknown field 
within each rod in a basis of cylindrical wave functional 
Owing to the inherent wave nature of the basis functions, 
each with the frequency (b guess , this basis is particularly 
well suited for the expansion, and the integral formula¬ 
tion in terms of the Green tensor means that the radi¬ 
ation condition is inherently satisfied. In total, this en¬ 
ables very high accuracy calculations of QNMs, although 
the method _is best suited to collections of cylindrical 23 
or spherical scatterers for which the projections onto 
the basis functions simplify significantly. 

To help in the following discussion, it is worth distin¬ 
guishing between the generalized effective mode volume, 
defined formally in Eq. © as a single complex number, 
and the numerically calculated numbers u" um (I?.), with 
i £ {Lai, Muljarov}, that is obtained from Eq. © by 
substituting the corresponding normalization formulas 
calculated in a circular domain of radius R. The up¬ 
per panel of Fig. [3] shows a zoom in at the real part 
of which an oscillatory behavior is clearly 

visible; from the discussion in Section |II A[ we can as¬ 
sociate these oscillations with the residual I?-dependent 
term A minimum in the oscillation amplitudes 

occurs at R ss 15a, after which the amplitudes increase 
without bounds as R —>• oo. 

An illustrative alternative display of the oscillatory be¬ 
havior of is found by plotting the values in the 




Numerical Integral Re{uj)™}/a 2 

FIG. 3: (Color online) Top: Real part of the calculated gen¬ 
eralized effective mode volume as function of calculation do¬ 
main radius R. Blue and red parts of the curve indicate the 
regions in which the oscillations decrease and increase in mag¬ 
nitude, respectively. Black dashed curves show the second 
order running averages starting from the first local maximum 
and minimum. Bottom: Values of the calculated generalized 
effective mode volume as a parameterized function of calcula¬ 
tion domain radius R for the part of the curve with decreasing 
oscillation amplitudes, 2 a < R < 15 a. The center of the in¬ 
nermost circle corresponds to the generalized effective mode 
volume and is indicated with a cross. 


complex frequency plane as shown in the lower panel 
of Fig. [3j In this way one can quite easily appreci¬ 
ate how the values of u£N n (I?) spirals initially towards 
a fixed value until the complex exponential factor in 
fhai(R) takes over and causes the values to spiral out¬ 
wards. The oscillating behavior provides a convenient 
and very graphical means of regularizing the integral by 
simply determining the center of the spiral. Clearly, this 
can be done with an accuracy which is much better than 
the smallest oscillation magnitude (better than the ra¬ 
dius of the innermost circle in the spiral). We find that 
the qualitative behavior of t>j™“(.R) in the form of a spi¬ 
ral is not affected by the accuracy of the calculations, 
but the center of the spiral is. Therefore, for all practical 
purposes, the error in the calculated value of vq, in this 
case, is dominated by the numerical integration and not 
the regularization of the integral. 

In practical calculations we find that one can get a 
convenient high-accuracy estimate of the center of the 
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spiral by calculating the moving average 


1 N 

V Q (Rn) = ^ V Q _1) ( jR «)> ( 31 ) 

n—riQ 

where Vq \R n ) denotes the ra’th calculated value of 
In the present case of the photonic crystal cav- 
ity, the second order moving average Vq shows a smooth 
behavior for 0 < R/a < 40, and by choosing two differ¬ 
ent starting points no, corresponding to consecutive local 
maxima and minima of the curve of nj™™ (/?), we get cor¬ 
responding upper and lower bounds for vq , as illustrated 
in Fig. [3j Because of the exponentially increasing factor 
in the running averages will in general not tend 

to a fixed value at large distances, but may eventually 
start increasing (exponentially). Therefore, we do not 
consider their use to be a formal regularization procedure 
in itself, but rather a simple computational method for 
estimating the center of the spiral. With this procedure 
we find the value 


vq/ a 2 = 0.988918 - 0.091688i, (32) 

with an estimated error of \5vq\/(i 2 < 2 x 10~ 6 . The os¬ 
cillating behavior in Fig. [3] is exactly what was noted in 
Ref. [5] which stated that “For very low-Q cavities, how¬ 
ever, the convergence is nontrivial due to the exponential 
divergence of the modes that may cause the inner prod¬ 
uct to oscillate around the proper value as a function of 
calculation domain size”. In practice, the magnitude of 
the oscillations are often much smaller than the desired 
accuracy for positions within typical calculation domains, 
and may also be much smaller than the overall numerical 
accuracy. 


Next, we compare to results of calculations based on 
the formulation by Muljarov et al.. Notably, in Eq. (18) 
there is no requirement of an infinite calculation do¬ 
main, which is very convenient from a numerical point 
of view. For the practical implementation, however, one 
must evaluate (numerically) both the first and the sec¬ 
ond derivative; together with the general accuracy of the 
integral, the accuracy of the derivatives then governs the 
overall precision. Because of the finite accuracy in the 
calculations, there is a numerically induced phase differ¬ 
ence between the two terms in Eq. (18), wherefore the 
changes in the volume integral is not fully compensated 
by the surface term. This difference becomes more pro¬ 
nounced at larger values of R , because the derivatives 
become larger. Instead of investigating the behavior of 
Eq. (181 as a function of calculation domain size as in 
^g. nr we take a fixed calculation domain size of R = 2a 
and evaluate the integrals with a fixed relative accuracy 
a 2 = 10 -7 and first order finite difference 


Muljarov I 


of 1 5' 

approximations to the derivatives calculated using suc¬ 
cessively smaller differences Ar. The results are shown 
in Table [I] and display a convergent behavior towards a 
value consistent with the calculation using the formula¬ 
tion of Lai et al ., cf. Eq. (32). 


A r/a ^iX? jarov (2 a)/a 2 

0.01 0.988894 - 0.091660i 

0.001 0.988918 - 0.091688i 
0.0001 0.988918 - 0.091688i 

TABLE I: Values of the calculated generalized effective mode 
volume u{X5) arov (2a) for a fixed calculation domain size R = 
2 a and relative integration tolerance <5r>M)ry a rov/ fl2 = 10 _ ' but 
varying accuracy of the numerical derivatives. 


Last, we compare to calculations based on the formu¬ 
lation by Sauvan et al. in Eq. (14). Our implemen¬ 
tation of the Fredholm integral equation does not al¬ 
low for evaluation of the QNM at complex positions. 
Instead, we used finite-difference time-domain (FDTD) 
calculationd- 25 ! with a run-time Fourier transform to cal¬ 
culate the QNM in the entire calculation domain, includ¬ 
ing the PML. The use of FDTD comes with a relatively 
large numerical error compared to the Fredholm inte¬ 
gral equation approach, but has the advantage that it 
is immediately applicable to completely general material 
structures. For these calculations, we used the formu¬ 
lation in Eq. (14) based on both electric and magnetic 
fields, and we evaluated the components of the fields 
at their individual positions in the Yee cells so as to 
avoid interpolation errors. In principle, the evaluation 
of Eq. (14), when calculated through the PML, should 
be independent of the calculation domain size. In prac¬ 
tice, however, we see a small variation with domain size, 
which we attribute to the limited accuracy of the numer¬ 
ical data. Table HI] summarizes the results calculated at 
four different domain sizes; in all cases the relative er¬ 
ror, when compared to the value in Eq. (32), is less than 
2%. Finally, using the same discretization, we can also 
use the FDTD data to calculate the generalized effective 
mode volume with the formulation of Lai et al. for which 
we find u™ TD /a 2 = 0.988 — 0.092i correspondning to a 
relative error of |<5uq|/|uq| < 0.001 when comparing to 
Eq. ([32]). 


L/a VsZLJa 2 

6 0.999 - 0.091i 

8 0.979 - 0.084i 

10 0.991 - 0.102i 

12 0.992 - 0.085i 

TABLE II: Values of the calculated generalized effective mode 
volume Ug™ an for different sizes of a square calculation do¬ 
main with side length L. 













C. Plasmonic nanorod dimer 


As a final example, we consider the three dimensional 
metallic nanorod dimer in Fig. |TJ which was recently 
studied in Ref. [22] It consists of two gold nanorods of 
length 100 nm and radius r ro d = 15 nrn in a homoge¬ 
neous background with refractive index tib = 1.5. The 
gold nanorods are aligned along the y axis and spaced by 
a gap of 20 nm. The relative permittivity of the gold is 
modeled as 

W D 

CrodM = 1- 2 ~T- -, (33) 

or + lwy 

with w p = 1.26 x 10 16 rad/s, and 7 = 1.41 x 10 14 rad/s. 
The resonance frequency of the dipolar like QNM of in¬ 
terest is uj c /2tt = (291 — 20i) THz, corresponding to the 
wavelength A c = 2'kc/lo = (1025 + 70i) nm and a quality 
factor of Q « 7. 

These QNM calculations were also performed using 
FDTE^with a run-time Fourier transform based on the 
raw Yee cell data without interpolation, which can be 
particularly detrimental in plasmonic structures due to 
the large held gradients. For the excitation, we used a 
spatial plane wave excitation near 291 THz with polar¬ 
ization along the dimer axis and a 6 fs Gaussian tempo¬ 
ral profile. As is common practise in FDTD when an¬ 
alyzing metals and other high index materials, we used 
a non-uniform mesh consisting of relatively small cubic 
Yee cells with side length 0.8 nm along the dimer axis and 
1 nm in the x and z directions for the inner region and 
larger cells of side length 50 nm for the outer region. See 
Ref. (223 for further calculation details. Despite the fine 
discretization, we expect the overall accuracy of these cal¬ 
culations to be limited to a few percent due to meshing, 
sub-meshing and numerical dispersion in FDTD, which 
is known to be particularly demanding for metallic struc¬ 
tured. We consider these to be state-of-the-art gridding 
parameters for metal resonators using FDTD, with sim¬ 
ulation times on the order of days on a small computer 
cluster. 

As in the case of the two-dimensional cavity, we can 
investigate the variation in the calculated generalized ef¬ 
fective mode volume with calculation domain size. To 
this end, we define the center position r c as the point 
directly between the two nanorods and evaluate Eq. |5]), 
with the substitution e r (r) —» cr(r, c5p) , in a rectangular 
cuboid of side lengths L x = L z and L y = 7 qm com¬ 
pletely enclosing the dimer. Figure [4] shows the variation 
in as a function of domain width L x . The varia¬ 

tion in the numerical integral shows a clear oscillation 
as in Fig. [2j although the spiraling behavior is less sym¬ 
metric. Despite the lower symmetry of the spiral, which 
we attribute to the lower degree of symmetry in the cal¬ 
culations, we can use the second order running averages 
to determine the center of the oscillation with relatively 
high precision. With this approach, we find 

' y Q/ r rod = 36.15 — 0.81i. (34) 




35.9 36 36.1 36.2 36.3 36.4 

Numerical Integral Re-fuJ^i 11 }/^ , 


FIG. 4: (Color online) Top: Real part of the calculated gener¬ 
alized effective mode volume as function of calculation domain 
width L x . Black dashed curves show the second order running 
averages. Bottom: Values of the calculated generalized effec¬ 
tive mode volume as a parameterized function of calculation 
domain width. The full curve is a spline-interpolation of the 
calculated data points shown by gray circles. The generalized 
effective mode volume vq is indicated with a cross. 


From the curves of the second order running averages, we 
estimate that the error stemming from the regularization 
procedure is |i5uq| < 0.01rj) od . This estimate does not 
take into account the error in the FDTD data or the 
error from the numerical quadrature. As noted above, 
we expect that the actual uncertainty in the FDTD data 
may be as large as a few percent, wherefore the accuracy 
in vq also in this three dimensional calculation is limited 
by the accuracy of the numerical calculations and not the 
regularization. 

Next, we compare to calculations based on the formu¬ 
lation by Sauvan et al. in Eq. (14). With this approach, 
we found vq/t 3 = 36.17— 1.02i, corresponding to a rela¬ 
tive difference of |Auq|/|cq| < 0.006 when comparing to 
Eq. (34) . In practise, we have found that the formulation 
in Eq. (14) is more sensitive to boundary effects in the 
FDTD data than Eq. ©• The difference is likely caused 
by the material weight function being larger inside the 
metal resonator, |cr(r,w)/??(r,u;)| <C 1, cf. Eqs. © and 
(14). In particular, we stress that one should take care 
not to use interpolated field values for the integration, 
but rather use the raw data at the different spatial grids 
for the individual components of the fields. 

Last, we remark on our attempts of using Eq. ( fl8| ) for 
calculating the QNM norm of the plasmonic dimer. We 
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have found the normalization of Muljarov et al. very 
difficult to evaluate with the FDTD data to a precision 
in which it provides better accuracy than Eq. ©• We 
believe that these problems are caused by the finite ac¬ 
curacy of the QNM calculation, which in practice is an 
unfortunate but rather general characteristic of many nu¬ 
merical Maxwell solvers in three dimensions. For FDTD, 
the inherent formulation in terms of Cartesian coordi¬ 
nates leads to additional interpolation errors when per¬ 
forming the radial differentiation as dictated by Eq. (181. 
In the present case of the plasmonic dimer, we found that 
the field gradients lead to relatively large numerical er¬ 
rors in the surface integral which in practice dominates 
the achievable accuracy. From the analysis in Section 
PH it is obvious that these limitations are not principal 
in nature, but rather practical issues pertaining to the 
numerical calculation method. 


Discussion 


Comparing Eqs. and ( p~8] ) , we can understand 

the three formulations of the norm as arising from three 
different approaches to the regularization of the integral 


j e r (r)f M (r)-f M (r)dV. 


(35) 


This integral can be regularized immediately by a com¬ 
plex coordinate transform, as discussed in section |II A| 
Therefore, in view of Eqs. (181 one may view the surface 
term in Eq. (J5| as a simple choice of regularization for 
use with data at positions along the real axis. It is worth 
noting, however, that the surface term appears dir ectl y 
in perturbation theory calculation^®, or from Eq. ( |17| , 
under the assumption of the Silver-Miiller radiation con¬ 
dition, and it is in principle important for the argument 
based on Eqs. (101 or (ll7|) that one can assign a well- 


defined value to the norm as R —> oo. In many practical 
calculations we find that the relatively simple formula 
of Lai et al. is very often sufficient, and that the un¬ 
certainty in the norm (or the associated effective mode 
volume) is dominated by the uncertainty in the numeri¬ 
cal data. When one has access to complex position data, 
the normalization of Sauvan et al. represents an obvi¬ 
ous alternative, and for high-accuracy numerical data, or 
when the QNMs are known analytically, one can benefit 
the most from the formula of Muljarov et al.. 

In this Article, we have focused on normalization of 
resonators embedded in a homogeneous background ma¬ 
terial, for which we have argued that Eq. ([2]) represents 


the proper choice of radiation condition. For resonators 
coupled to waveguides the radiation condition is different, 
but also in this case there appears to be an inherent need 
for regularization of the normalization integral 1 2 '. Also, 
it is worth noting that Bai et al. recently introduced 
an alternative calculation method for QNMs, based on a 
scattering formulation at complex frequencies, in which 
the resulting QNMs are naturally normalized without an 
explicit need for an integration stepP^. 


IV. CONCLUSIONS 


In conclusion, we have discussed the relation between 
three different formulations for the normalization of 
QNMs in leaky optical cavities and plasmonic res¬ 
onators. The three formulations arise from independent 
derivations of an expansion for the electromagnetic 
field in terms of QNMs, and for spherical calculation 
volumes, we have shown the equivalence between 
the different formulas explicitly. Moreover, we have 
illustrated the equivalence between the norms with a 
number of example calculations of practical interest. 
We have argued, that the apparent differences in the 
normalization formulas can be understood as arising 
from alternative ways of regularizing the normalization 
integral to handle the inherently divergent QNMs. In 
this view, the surface term in the normalization by 
Lai et al. represents a simple choice of regularization 
which in principle is insufficient to properly regularize 
the integral when the size of the calculation domain is 
varied along the real axis. We have discussed this issue 
in detail, and we have shown how one can, in principle, 
always regularize the integral by a complex coordinate 
transform. In practice, we find that there is rarely any 
need for additional regularization beyond the relatively 
simple formula of Lai et al. Nevertheless, depending 
on the calculation method used to obtain the QNMs, 
the normalizations of Sauvan et al. or Muljarov et al. 
may be more convenient. Regardless of the choice of 
normalization integral, we have discussed how one can 
use the norm to define an effective mode volume for use 
in the Purcell formula. Because of the complementarity 
of the three norms, this naturally leads to the same 
estimates of the enhanced spontaneous emission factor. 
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